// standard headers
#include <istream>
#include <vector>
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <map>
#include <cstdlib>

// boost headers
#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/filesystem/path.hpp>

// components headers
#include <print_error.h>
#include <tree_kitsch.h>
#include <tree.h>
#include <tree_ops.h>
#include <portable_timer.h>
#include <progress_indicator.h>
#include <matrix.h>
#include <matrix_io.h>
#include "get_arg.h"
#define DEBUG_LEVEL 6
#include <debug_func.h>




using boost::format;
using boost::lexical_cast;
using boost::filesystem::path;
using std::vector;
using std::ostringstream;

using std::string;
using std::fill;
using std::ofstream;
using std::cout;
using std::cin;
using std::cerr;
using std::ifstream;
using std::map;
using std::pair;
using std::istream;
using std::max;
using std::min;
using namespace bcpp_matrix;

struct t_incorrect_pairs
{
	string cluster_name;
	unsigned cnt_genes;
	unsigned cnt_missing_pairs;
	unsigned cnt_expected_pairs;
	t_incorrect_pairs(const t_matrix_f& data, const string& cluster_name_):
		cluster_name(cluster_name_),
		cnt_genes(data.size()),
		cnt_missing_pairs(data.get_count_missing()),
		cnt_expected_pairs(data.get_count_pairs())
	{}
};


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//	main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

int main (int argc, char *argv[])
{
	//	we use cout to print progress so flush per call
	//	cout.setf(std::ios_base::unitbuf);

	std::ios::sync_with_stdio(false);

	t_program_args args;
	try
	{
	if (!process_command_line_options(argc, argv, args))
		return 1;



	while (cin)
	{
		//
		//	save clusters with missing entries
		//
		vector<t_incorrect_pairs> problems;
	
	
	
		//
		//	read data from input file
		//
		vector<t_matrix_f>	data;
		args.VERBOSE() << "\tKitsch\n";
		args.VERBOSE() << "\tReading matrix and single linkage clustering...\n";
		unsigned cnt_lines = 0;
		string tree_name;
		{
	
			t_progress_indicator dots(args.VERBOSE(), 20000);
			t_matrix_f initial_matrix;

			cin >> t_matrix_istrm_wrapper_f(initial_matrix, args.input_format, cnt_lines, dots);
			if (!cin)
				return 0;
			//if (!bcpp_matrix::read_matrix(cin,  initial_matrix,  args.input_format, cnt_duplicate_pairs, dots);
			tree_name = initial_matrix.matrix_name;
			if (!tree_name.length())
				tree_name = "kitsch optimized tree";
			initial_matrix.get_single_linkage_clusters(data);
		}
	
	
		if (!data.size())
		{
			std_print_warning("Empty matrix supplied.");
			return 1;
		}
	
	
	
		//
		//	iterate through each subcluster
		//
		std::vector<string> tree_strings;
		for (unsigned i = 0; i < data.size(); ++i)
		{
			string cluster_id = (data.size() > 1) ? "cluster " + lexical_cast<string>(i + 1) : "";
			if (!data[i].matrix_name.length() && data.size() > 1)
				data[i].matrix_name = cluster_id;
			t_matrix_f& sub_cluster = data[i];
			// save problems
			if (sub_cluster.get_count_missing())
				problems.push_back(t_incorrect_pairs(sub_cluster, data[i].matrix_name));
	
			args.VERBOSE()	<< "\t" << tree_name
							<< cluster_id <<" has " << data[i].size() << " elements\n";
	
	
			bcpp_tree::t_tree kitsch_tree;
			//888888888888888888888888888888888888888888888888888888888888888888888888
	
			// kitsch
	
			//888888888888888888888888888888888888888888888888888888888888888888888888
			{
				args.VERBOSE() << "\tBuilding kitsch tree...\n";
				t_progress_indicator dots(args.VERBOSE(), 20000);
				dots.use_timer();
				unsigned							best_jumble;
				float								sum_sqrs;
				float								avg_percent_stddev;
				unsigned							cnt_trees_examined;
				kitsch_tree::kitsch_make_tree(kitsch_tree, sub_cluster.names, sub_cluster.matrix,
												args.cnt_iterative_fitch_passes,
												args.power,
												best_jumble,
												sum_sqrs,
												avg_percent_stddev,
												cnt_trees_examined, dots);
				args.VERBOSE() << format("\t%|1$4| trees examined\n"
										"\tBest results from iteration %|2$| out of %|3$|\n"
										"\tSum of squares         = %|4$.5|\n"
										"\tAverage percent stddev = %|5$.5f|\n")
						% cnt_trees_examined
						% best_jumble
						% args.cnt_iterative_fitch_passes
						% sum_sqrs
						% avg_percent_stddev;
	
			}
	
	
			ostringstream ostrm;
			bcpp_tree_ops::do_newhampshire_format_str nh_str(ostrm);
			kitsch_tree.for_each(nh_str);
			tree_strings.push_back(ostrm.str());
			// remove last ";"
			tree_strings.back().erase(tree_strings.back().length() - 1);
		}
	
		// get tree in new_hampshire format
		if (args.header)
			cout << ">" << tree_name << "\n";
		if (tree_strings.size() == 1)
			cout << tree_strings[0] << ";\n";
		else
		{
			args.VERBOSE() << "\tJoining " << tree_strings.size() << " unrelated sub-trees...\n";
			cout << "(";
			cout << tree_strings[0] << ":0.0";
			for (unsigned i = 1; i < tree_strings.size(); ++i)
				cout << "," << tree_strings[i] << ":0.0";
			cout << ");\n";
		}
	
	
		//
		//	print clusters with missing entries
		//
		if (problems.size())
		{
			std_print_warning(args.VERBOSE(),  "Some matrixes had missing pairs.");
			for (unsigned i = 0; i < problems.size(); ++i)
			{
	
				args.CERR() << format("\t\t%|1$4|: "
										"%|2$4| genes, "
										"%|3$4| gene pairs expected, "
										"%|4$4| gene pairs missing.\n ")
						% problems[i].cluster_name
						% problems[i].cnt_genes
						% problems[i].cnt_expected_pairs
						% problems[i].cnt_missing_pairs;
			}
		}
		args.VERBOSE() << "\tAll finished!!\n\n";
	}
	}
	catch(std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

